#ifndef __ARRAY2D_H__
  #define __ARRAY2D_H__

#include <string>
#include <sstream>

#include "exceptions.h"
#include "macros.h"
#include "global_variables.h"

using namespace std;

//! Two dimensional array

/*!
A simple array implementation.
For arrays larger than 1GB some changes may be necessary.

Implementation date: 2008

Usage:

\code
array2d<int> a(2, 3, 0); // with 2 rows and 3 columns all initialized with 0

a(0, 0) = 17;

a(0, 1) = 21;

a(1, 0) = -4;
\endcode
*/

template <class T>
class array2d {
private:
  T *content;
  unsigned int stream_position;
  unsigned int rows;
  unsigned int cols;

  void allocate();
  void release();
  void hard_copy(const array2d &a);

public:
  //! resets the array to an empty array and releases allocated memory
  void clear();

  //! initializes all elements of the array to init_value
  void initialize(const T &init_value);

  //! returns the position in which the next element will be inserted when operator<< is called
  unsigned int get_stream_position() const;

  //! sets the position in which the next element will be inserted when operator<< is called
  void set_stream_position(unsigned int position);

  //! returns the string representation of the object
  string to_string() const;

  //! returns the number of rows
  unsigned int get_row_count() const;

  //! returns the number of columns
  unsigned int get_column_count() const;

  //! reduces the size of the array (may throw an exception)
  void cut_rows_after(unsigned int index);

  //! reconstructs the array and deletes its content
  void reconstruct(unsigned int rows, unsigned int cols);

  //! resizes the array and keeps its content intact
  void resize(unsigned int rows, unsigned int cols, const T &init_value);

  //! negates and returns a copy of the array
  array2d negate() const;

  //! mathematical addition (may throw an exception)
  array2d operator+(const array2d &a) const;

  //! mathematical subtraction (may throw an exception)
  array2d operator-(const array2d &a) const;

  //! mathematical multiplication
  array2d operator*(const T &value) const;

  //! mathematical multiplication
  array2d operator*(const array2d &a) const;

  //! index operator (may throw an exception)
  T &operator()(unsigned int row, unsigned int col) const;

  //! assignment operator
  array2d &operator=(const array2d &a);

  //! shift operator: inserts an element into the array at the stream_position (may throw an exception)
  array2d<T> &operator<<(const T &value);

  //! copy constructor
  array2d(const array2d &a);

  //! default constructor: an empty array, reconstruct or resize must be called to enlarge the array before it can be used
  array2d();

  //! constructor: does not initialize the array content
  explicit array2d(unsigned int rows, unsigned int cols);

  //! constructor: initializes the array content
  explicit array2d(unsigned int rows, unsigned int cols, const T &init_value);

  //! destructor (non-virtual: please do not inherit from this class)
  ~array2d();
}; // class array2d

//////////////////////////////////////////////////////////////////////////////

template <class T>
void array2d<T>::allocate() {
  release();
  content = new T[rows * cols];
}

template <class T>
void array2d<T>::release() {
  delete[] content;
  content = 0;
}

template <class T>
void array2d<T>::hard_copy(const array2d &a) {
  this->stream_position = 0;
  this->rows = a.rows;
  this->cols = a.cols;

  allocate();
  const unsigned int total = rows * cols;

  #pragma omp parallel for if (total > MINIMUM_SIZE_LIMIT_OMP)
  for (unsigned int i = 0; i < total; ++i) // do not use memcpy (it may cause problems with non-standard datatypes)
    this->content[i] = a.content[i];
}

template <class T>
void array2d<T>::clear() {
  release();
  stream_position = 0;
  rows = 0;
  cols = 0;
}

template <class T>
void array2d<T>::initialize(const T &init_value) {
  const unsigned int total = rows * cols;

  #pragma omp parallel for if (total > MINIMUM_SIZE_LIMIT_OMP)
  for (unsigned int i = 0 ; i < total; ++i)
    content[i] = init_value;
}

template <class T>
unsigned int array2d<T>::get_stream_position() const {
  return this->stream_position;
}

template <class T>
void array2d<T>::set_stream_position(unsigned int position) {
  this->stream_position = position;
}

template <class T>
string array2d<T>::to_string() const {
  ostringstream result;

  for (unsigned int i = 0 ; i < rows; ++i) {
    for (unsigned int j = 0 ; j < cols; ++j) {
      result << '[';
      result << i;
      result << ", ";
      result << j;
      result << "] = ";
      result << content[i * this->cols + j];
      result << '\n';
    } // for j
  } // for i

  return result.str();
}

template <class T>
unsigned int array2d<T>::get_row_count() const {
  return rows;
}

template <class T>
unsigned int array2d<T>::get_column_count() const {
  return cols;
}

template <class T>
void array2d<T>::cut_rows_after(unsigned int index) {
  rows = min(index, rows);
}

template <class T>
void array2d<T>::reconstruct(unsigned int rows, unsigned int cols) {
  this->rows = rows;
  this->cols = cols;

  allocate();

  set_stream_position(0); // makes sure that operator<< will begin from the index 0
}

template <class T>
void array2d<T>::resize(unsigned int rows, unsigned int cols, const T &init_value) {
  if (this->rows >= rows && this->cols >= cols)
    return; // do not shrink the array now!

  const unsigned int total = rows * cols;
  T *tmp = new T[total];

  #pragma omp parallel for if (rows > MINIMUM_SIZE_LIMIT_OMP)
  for (unsigned int i = 0 ; i < rows; ++i) {
    for (unsigned int j = 0 ; j < cols; ++j) {
      if (i < this->rows && j < this->cols)
        tmp[i * cols + j] = content[i * this->cols + j];
      else
        tmp[i * cols + j] = init_value;
    } // for j
  } // for i

  this->rows = rows;
  this->cols = cols;

  release();
  content = tmp;
}

template <class T>
array2d<T> array2d<T>::negate() const {
  return (*this) * -1;
}

template <class T>
array2d<T> array2d<T>::operator+(const array2d &a) const {
  if (a.rows != this->rows || a.cols != this->cols)
    throw operation_not_allowed_exception(global_empty_string);

  const unsigned int total = rows * cols;
  array2d result(*this);

  #pragma omp parallel for if (total > MINIMUM_SIZE_LIMIT_OMP)
  for (unsigned int i = 0 ; i < total; ++i)
    result.content[i] += a.content[i];

  return result;
}

template <class T>
array2d<T> array2d<T>::operator-(const array2d &a) const {
  return (*this + a.negate());
}

template <class T>
array2d<T> array2d<T>::operator*(const T &value) const {
  array2d result(*this);
  const unsigned int total = rows * cols;

  #pragma omp parallel for if (total > MINIMUM_SIZE_LIMIT_OMP)
  for (unsigned int i = 0 ; i < total; ++i)
    result.content[i] *= value;

  return result;
}

template <class T>
array2d<T> array2d<T>::operator*(const array2d &a) const {
  if (this->cols != a.rows)
    throw operation_not_allowed_exception(global_empty_string);

  array2d result(this->rows, a.cols);

  #pragma omp parallel for if (rows > MINIMUM_SIZE_LIMIT_OMP)
  for (unsigned int i = 0; i < this->rows; ++i) {

    T tmp;
    for (unsigned int j = 0; j < a.cols; ++j) {
      tmp = 0;

      for (unsigned int k = 0; k < this->cols; ++k)
        tmp += this->operator()(i, k) * a(k, j);

      result(i, j) = tmp;
    } // for j

  } // for i

  return result;
}

template <class T>
T &array2d<T>::operator()(unsigned int row, unsigned int col) const {
  INDEX_CHECK(row, this->rows);
  INDEX_CHECK(col, this->cols);

  return content[row * this->cols + col];
}

template <class T>
array2d<T> &array2d<T>::operator=(const array2d &a) {
  if (&a == this)
    return *this;

  hard_copy(a); // content will be released first
  return *this;
}

template <class T>
array2d<T> &array2d<T>::operator<<(const T &value) {
  INDEX_CHECK_MSG(stream_position, rows * cols, "Out of range! enlarge the array first.");

  content[stream_position] = value;
  ++stream_position;

  return *this;
}

template <class T>
array2d<T>::array2d(const array2d &a) : content(0), stream_position(0), rows(a.rows), cols(a.cols) {
  hard_copy(a);
}

template <class T>
array2d<T>::array2d() : content(0), stream_position(0), rows(0), cols(0) {
  // nothing to do
}

template <class T>
array2d<T>::array2d(unsigned int rows, unsigned int cols) : content(0), stream_position(0), rows(rows), cols(cols) {
  allocate();
}

template <class T>
array2d<T>::array2d(unsigned int rows, unsigned int cols, const T &init_value) : content(0), stream_position(0), rows(rows), cols(cols) {
  allocate();
  initialize(init_value);
}

template <class T>
array2d<T>::~array2d() {
  release();
}

#endif
